# Replication File for
# "Coalition Prospects and Policy Change: An Application to the Environment"
# Authors: Mark A. Kayser and Jochen Rehmert


# load required packages
library(ggplot2);library(grid);library(stargazer);library(plm);library(lmtest);library(sandwich)

# set directory
setwd('C:/Users/rehmerjo/Dropbox/CIP_LSQ_Conference/dataverse')

# read data
dat <- read.csv("EPS_data.csv", sep = ";", dec = ",")

# NB: The file has been saved on a German machine.
# This means that any .csv file will use the ";" as 
# a delimiter and the "," as the decimal separator

# coding of control variables 
dat$kyoto <- 0
dat$kyoto[dat$country == "austria" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "czechia" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "denmark" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "estonia" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "finland" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "germany" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "greece" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "ireland" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "italy" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "japan" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "netherlands" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "new zeland" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "norway" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "poland" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "portugal" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "slovakia" & dat$year >= 1999] <- 1
dat$kyoto[dat$country == "slovenia" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "spain" & dat$year >= 1998] <- 1
dat$kyoto[dat$country == "sweden" & dat$year >= 1998] <- 1

dat$fukushima <- 0
dat$fukushima[dat$year >= 2011] <- 1

##################################################################
# Figure 1                                                       #
##################################################################
# read in additional data on monthly polling/CIP 
dat2 <- read.csv("Germany_MerkelII.csv", sep = ";", dec = ",")

dat2$date <- as.Date(dat2$date, "%d.%m.%Y")

dat2 <- dat2[order(dat2$date),]

par(mfrow = c(2,1),
    oma = c(2,3,0,0)+0.1,
    mar = c(0,0,1,1) + 0.1,
    mai = c(0.1,0.1, 0.1,0.1),
    mgp = c(2,1,0))
plot(dat2$cip_gross[dat2$party_abbr == "CDU-CSU"] ~ as.Date(dat2$date[dat2$party_abbr == "CDU-CSU"], origin = "1899-12-31")
     , ylim = c(0,.6), ylab = "CIP (net)", xaxt='n', col = "white", xlab = "")

lines(dat2$polls[dat2$party_abbr == "FDP"]/100 ~ as.Date(dat2$date[dat2$party_abbr == "FDP"], origin = "1899-12-31"), col=colors()[654], lty = 1, lwd = 4)
abline(h = 0.05, lty = 2)
text(x = as.Date("2010-05-01"), y = 0.025, label = "5% Hurdle   FDP's Polls")

lines(dat2$polls[dat2$party_abbr == "CDU-CSU"]/100 ~ as.Date(dat2$date[dat2$party_abbr == "CDU-CSU"], origin = "1899-12-31"), lty = 1, lwd = 3)
text(x = as.Date("2012-06-15"), y = 0.45, label = "CDU/CSU's Polls")

lines(dat2$polls[dat2$party_abbr == "Greens"]/100 ~ as.Date(dat2$date[dat2$party_abbr == "Greens"], origin = "1899-12-31"), col=c("green3"),lty = 1, lwd = 3)
text(x = as.Date("2012-01-01"), y = 0.2, label = "Greens' Polls")

mtext("Poll Average", side = 2, line = 2, cex = 1.5)


plot(dat2$cip_eco[dat2$party_abbr == "CDU-CSU"] ~ as.Date(dat2$date[dat2$party_abbr == "CDU-CSU"], origin = "1899-12-31")
     , ylim = c(.4,1), ylab = "CIP (net)", xaxt='n', col = "white", xlab = "")
lines(dat2$cip_eco[dat2$party_abbr == "CDU-CSU"] ~ as.Date(dat2$date[dat2$party_abbr == "CDU-CSU"], origin = "1899-12-31"), lty = 1, lwd = 2)
segments(as.Date(dat2$date[dat2$party_abbr == "CDU-CSU"], origin = "1899-12-31"), as.numeric(dat2$cip_eco_lower[dat2$party_abbr == "CDU-CSU"]), as.Date(dat2$date[dat2$party_abbr == "CDU-CSU"], origin = "1899-12-31"), as.numeric(dat2$cip_eco_upper[dat2$party_abbr == "CDU-CSU"]))
segments(as.Date(dat2$date[dat2$party_abbr == "CDU-CSU"]-2, origin = "1899-12-31"), as.numeric(dat2$cip_eco_lower[dat2$party_abbr == "CDU-CSU"]), as.Date(dat2$date[dat2$party_abbr == "CDU-CSU"] + 2, origin = "1899-12-31"), as.numeric(dat2$cip_eco_lower[dat2$party_abbr == "CDU-CSU"]))
segments(as.Date(dat2$date[dat2$party_abbr == "CDU-CSU"]-2, origin = "1899-12-31"), as.numeric(dat2$cip_eco_upper[dat2$party_abbr == "CDU-CSU"]), as.Date(dat2$date[dat2$party_abbr == "CDU-CSU"] +2 , origin = "1899-12-31"), as.numeric(dat2$cip_eco_upper[dat2$party_abbr == "CDU-CSU"]))

text(x = as.Date(38000, origin = "1899-12-31" ), y = .65, label = "FDP")
text(x = as.Date(38000, origin = "1899-12-31" ), y = .8, label = "CDU/CSU")
text(x = as.Date(38000, origin = "1899-12-31" ), y = .235, label = "B90/Gr�ne")
text(x = as.Date(38000, origin = "1899-12-31" ), y = .425, label = "SPD")
text(x = as.Date(38000, origin = "1899-12-31" ), y = .05, label = "Die Linke")
mtext("CIP", side = 2, line = 2, cex = 1.5)
abline(v =as.Date("2011-03-11"), lty=1, lwd=3)
abline(v =as.Date("2011-05-30"), lty=1, lwd=3)
text(x = as.Date("2010-12-01"), y = .4, label = "Fukushima")
text(x = as.Date("2011-09-15"), y = .4, label = "Nuclear Phase-Out")
text(x = as.Date("2012-06-15" ), y = .93, label = "CDU/CSU's CIP excl. the Greens")

##################################################################
# Figure 2                                                       #
##################################################################

par(mfrow = c(3,3), mai = c(0.4, 0.4, 0.5, 0.1))
plot(NA, xlim = c(1990,2012), ylim = c(0,4), xlab= "", ylab = "", main = "Austria", cex.axis=1.5, cex.main = 1.5)
points(dat$eps[dat$country == "austria"] ~ dat$year[dat$country == "austria"])
lines(dat$eps[dat$country == "austria"] ~ dat$year[dat$country == "austria"])
plot(NA, xlim = c(1990,2012), ylim = c(0,4), xlab= "", ylab = "", main = "Denmark", cex.axis=1.5, cex.main = 1.5)
points(dat$eps[dat$country == "denmark"] ~ dat$year[dat$country == "denmark"])
lines(dat$eps[dat$country == "denmark"] ~ dat$year[dat$country == "denmark"])
plot(NA, xlim = c(1990,2012), ylim = c(0,4), xlab= "", ylab = "", main = "Finland", cex.axis=1.5, cex.main = 1.5)
points(dat$eps[dat$country == "finland"] ~ dat$year[dat$country == "finland"])
lines(dat$eps[dat$country == "finland"] ~ dat$year[dat$country == "finland"])
plot(NA, xlim = c(1990,2012), ylim = c(0,4), xlab= "", ylab = "", main = "Germany", cex.axis=1.5, cex.main = 1.5)
points(dat$eps[dat$country == "germany"] ~ dat$year[dat$country == "germany"])
lines(dat$eps[dat$country == "germany"] ~ dat$year[dat$country == "germany"])
plot(NA, xlim = c(1990,2012), ylim = c(0,4), xlab= "", ylab = "", main = "Ireland", cex.axis=1.5, cex.main = 1.5)
points(dat$eps[dat$country == "ireland"] ~ dat$year[dat$country == "ireland"])
lines(dat$eps[dat$country == "ireland"] ~ dat$year[dat$country == "ireland"])
plot(NA, xlim = c(1990,2012), ylim = c(0,4), xlab= "", ylab = "", main = "Italy", cex.axis=1.5, cex.main = 1.5)
points(dat$eps[dat$country == "italy"] ~ dat$year[dat$country == "italy"])
lines(dat$eps[dat$country == "italy"] ~ dat$year[dat$country == "italy"])
plot(NA, xlim = c(1990,2012), ylim = c(0,4), xlab= "", ylab = "", main = "Netherlands", cex.axis=1.5, cex.main = 1.5)
points(dat$eps[dat$country == "netherlands"] ~ dat$year[dat$country == "netherlands"])
lines(dat$eps[dat$country == "netherlands"] ~ dat$year[dat$country == "netherlands"])
plot(NA, xlim = c(1990,2012), ylim = c(0,4), xlab= "", ylab = "", main = "Portugal", cex.axis=1.5, cex.main = 1.5)
points(dat$eps[dat$country == "portugal"] ~ dat$year[dat$country == "portugal"])
lines(dat$eps[dat$country == "portugal"] ~ dat$year[dat$country == "portugal"])
plot(NA, xlim = c(1990,2012), ylim = c(0,4), xlab= "", ylab = "", main = "Sweden", cex.axis=1.5, cex.main = 1.5)
points(dat$eps[dat$country == "sweden"] ~ dat$year[dat$country == "sweden"])
lines(dat$eps[dat$country == "sweden"] ~ dat$year[dat$country == "sweden"])


##################################################################
# Table 1                                                        #
##################################################################

# run base-line model first to acquire row index of constant sample
summary(mod.1.a <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                       eco_cip_gross + eco_polls , 
                     data = dat[!is.na(dat$eco_polls) & dat$cee == 0 ,]))

# output descriptive statistics for those observations used in analyses
stargazer(dat[row.names(mod.1.a[["model"]]) ,c("eps","cab_minority", "kyoto","gdp_quarter","environ_cab_mean",
                 "total_ghg","eco_gov", "environ_eco","eco_cip_gross","pm_cip_eco","eco_polls")],
          covariate.labels = c("Environmental Policy Stringency","Minority Cabinet",
                               "Kyoto Protocol", "Quarterly GDP Growth (yearly mean)",
                               "Cabinet's Mean Environmental Protection","Total Greenhouse Gas Emissions/Capita",
                               "Green Party in Government","Green Party's Environmental Protection",
                               "Green Party's gross CIP (yearly mean)","PM Party's net CIP excl. Green Part(ies) (yr. mean)",
                               "Green Party's Polls (yearly mean)"), nobs = TRUE)


##################################################################
# Table 2                                                        #
##################################################################

### CIP ###
# fixed-effects, clustered SE
summary(mod.1.a <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                       eco_cip_gross , 
                     data = dat[!is.na(dat$eco_polls) & dat$cee == 0 ,]))
mod.1.a.se <- coeftest(mod.1.a, vcov=vcovHC(mod.1.a,type="HC0",cluster="group"))[,2]

# pooled with LDV
summary(mod.1.c <-lm(eps ~ l_eps + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                       eco_cip_gross + year , 
                     data = dat[!is.na(dat$eco_polls) & dat$cee == 0,]))
mod.1.c.se <- coeftest(mod.1.c, vcov=vcovHC(mod.1.c,type="HC0",cluster="group"))[,2]


# first-difference
summary(mod.1.b <-plm(eps ~  cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                        eco_cip_gross  , 
                      data = dat[!is.na(dat$eco_polls) & dat$cee == 0,], model = "fd", index = c("country_id","year")))
mod.1.b.se <- coeftest(mod.1.b, vcov=vcovHC(mod.1.b,type="HC0",cluster="group"))[,2]


### Polls ###
# fixed-effects, clustered SE
summary(mod.2.a <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                       I(eco_polls/100) , 
                     data = dat[  dat$cee == 0 & !is.na(dat$eco_cip_gross),]))
mod.2.a.se <- coeftest(mod.2.a, vcov=vcovHC(mod.2.a,type="HC0",cluster="group"))[,2]

# pooled with LDV
summary(mod.2.c <-lm(eps ~ l_eps + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                       I(eco_polls/100) + year , 
                     data = dat[!is.na(dat$eco_polls) & dat$cee == 0 & !is.na(dat$eco_cip_gross),]))
mod.2.c.se <- coeftest(mod.2.c, vcov=vcovHC(mod.2.c,type="HC0",cluster="group"))[,2]


# first-difference
summary(mod.2.b <-plm(eps ~  cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                        I(eco_polls/100)  , 
                      data = dat[!is.na(dat$eco_polls) & !is.na(dat$eco_cip_gross) & dat$cee == 0,], model = "fd", index = c("country_id","year")))
mod.2.b.se <- coeftest(mod.2.b, vcov=vcovHC(mod.2.b,type="HC0",cluster="group"))[,2]



stargazer(list(mod.1.a,   mod.1.c, mod.1.b, mod.2.a,   mod.2.c, mod.2.b),
          se = list(mod.1.a.se, mod.1.c.se, mod.1.b.se, mod.2.a.se, mod.2.c.se, mod.2.b.se),
          omit = c("country_id","year"),
          covariate.labels = c("EPS$_{t-1}$", "Minority Cabinet","Kyoto Protocol",
                               "GDP Growth","Cabinet's Mean Environmental Protection",
                               "Total Greenhouse Gas Emissions/Capita","Green Party in Government",
                               "Green Party's Environmental Protection",
                               "Green Party's gross CIP", "Green Party's Polls", "Intercept"),
          add.lines = list(c("Number of Countries","9","9","9","9","9", "9"),
                           c("Fixed-Effects", "Yes", "No","No", "Yes", "No","No"),
                           c("Lagged Dependent Vairable","No","Yes","No","No","Yes","No"),
                           c("Trend", "No", "Yes", "No","No", "Yes", "No") ), dep.var.labels = "")

##################################################################
# Table 3                                                        #
##################################################################

### CIP ###
# fixed-effects, clustered SE
summary(mod.3.a <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                       pm_cip_eco, 
                     data = dat[!is.na(dat$eco_polls) & dat$cee == 0,]))
mod.3.a.se <- coeftest(mod.3.a, vcov=vcovHC(mod.3.a,type="HC0",cluster="group"))[,2]


# pooled with LDV
summary(mod.3.c <-lm(eps ~ l_eps + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                       pm_cip_eco  + year , 
                     data = dat[!is.na(dat$eco_polls) & dat$cee == 0,]))
mod.3.c.se <- coeftest(mod.3.c, vcov=vcovHC(mod.3.c,type="HC0",cluster="group"))[,2]

# first-difference
summary(mod.3.b <-plm(eps ~  cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                        pm_cip_eco , 
                      data = dat[!is.na(dat$eco_polls) & dat$cee == 0,], model = "fd", index = c("country_id","year")))
mod.3.b.se <- coeftest(mod.3.b, vcov=vcovHC(mod.3.b,type="HC0",cluster="group"))[,2]

# fixed-effects, clustered SE & interaction 
summary(mod.3.d <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                       pm_cip_eco * eco_cip_gross, 
                     data = dat[!is.na(dat$eco_polls) & dat$cee == 0,]))
mod.3.d.se <- coeftest(mod.3.d, vcov=vcovHC(mod.3.d,type="HC0",cluster="group"))[,2]

# fixed-effects, clustered SE & interaction 
summary(mod.3.e <-lm(eps ~  as.factor(country_id) + cab_minority + kyoto+ gdp_quarter + environ_cab_mean + total_ghg + eco_gov +  environ_eco +
                       pm_cip_eco * eco_cip_gross + year, 
                     data = dat[!is.na(dat$eco_polls) & dat$cee == 0,]))
mod.3.e.se <- coeftest(mod.3.e, vcov=vcovHC(mod.3.e,type="HC0",cluster="group"))[,2]


stargazer(list(mod.3.a, mod.3.c, mod.3.b, mod.3.d, mod.3.e), 
          se = list(mod.3.a.se, mod.3.c.se, mod.3.b.se, mod.3.d.se, mod.3.e.se),
          omit = c("country_id","year"),
          covariate.labels = c("EPS$_{t-1}$", "Minority Cabinet","Kyoto Protocol",
                               "GDP Growth","Cabinet's Mean Environmental Protection",
                               "Total Greenhouse Gas Emissions/Capita","Green Party in Government",
                               "Green Party's Environmental Protection",
                               "PM Party's CIP net of Green parties (PM CIP)", "Green Party's gross CIP (Green CIP)", "PM CIP $\\times$ Green CIP", "Intercept"),
          add.lines = list(c("Number of Countries","9","9","9","9","9"),
                           c("Fixed-Effects", "Yes", "No","No","Yes","Yes"),
                           c("Lagged Dependent Vairable","No","Yes","No","No","No"),
                           c("Trend", "No", "Yes", "No", "No", "Yes") ), dep.var.labels = "")


##################################################################
# Figure 3, panel A first, then panel B                          #
##################################################################


plotDat1 <- data.frame(rbind(
  cbind(model = "LDV", estimate = mod.1.c$coefficients["eco_cip_gross"], se = mod.1.c.se["eco_cip_gross"], var = "CIP (gross)"),
  cbind(model = "LDV", estimate = mod.2.c$coefficients["I(eco_polls/100)"], se = mod.2.c.se["I(eco_polls/100)"], var = "Polls"   ),
  cbind(model = "Fixed-Effects", estimate = mod.1.a$coefficients["eco_cip_gross"], se = mod.1.a.se["eco_cip_gross"], var = "CIP (gross)"),
  cbind(model = "Fixed-Effects", estimate = mod.2.a$coefficients["I(eco_polls/100)"], se = mod.2.a.se["I(eco_polls/100)"], var = "Polls"   ),
  cbind(model = "First-Diff", estimate = mod.1.b$coefficients["eco_cip_gross"], se = mod.1.b.se["eco_cip_gross"], var = "CIP (gross)"),
  cbind(model = "First-Diff", estimate = mod.2.b$coefficients["I(eco_polls/100)"], se = mod.2.b.se["I(eco_polls/100)"], var = "Polls"   )))



plotDat1$estimate <- as.numeric(as.character(plotDat1$estimate))
plotDat1$se <- as.numeric(as.character(plotDat1$se))
plotDat1$upper95 <- plotDat1$estimate + plotDat1$se * qnorm(0.975)
plotDat1$lower95 <- plotDat1$estimate - plotDat1$se * qnorm(0.975)
plotDat1$upper90 <- plotDat1$estimate + plotDat1$se * qnorm(0.95)
plotDat1$lower90 <- plotDat1$estimate - plotDat1$se * qnorm(0.95)
plotDat1$var <- as.character(plotDat1$var)
plotDat1$model <- as.character(plotDat1$model)

dev.off()

par(mar = c(4,2,4,2) + 0.1, mai = c(1,1,.5,.5), oma = c(1,1,1,1))
plot(plotDat1$estimate[plotDat1$var == "CIP (gross)"] ~ c(1.5,1,2), ylim = c(0, 4), pch = 19, cex = 1.25, col = "black", xaxt='n', yaxt = 'n',  xlab = "",
     ylab = "Marginal Effect on EPS" , cex.lab = 1.7, frame = FALSE, xlim = c(0.95,2))
segments(1, plotDat1$estimate[plotDat1$var == "CIP (gross)" & plotDat1$model == "Fixed-Effects"], 1, plotDat1$lower95[plotDat1$var == "CIP (gross)" & plotDat1$model == "Fixed-Effects"], lty = 2)
segments(1, plotDat1$estimate[plotDat1$var == "CIP (gross)" & plotDat1$model == "Fixed-Effects"], 1, plotDat1$upper95[plotDat1$var == "CIP (gross)" & plotDat1$model == "Fixed-Effects"], lty = 2)
segments(2, plotDat1$estimate[plotDat1$var == "CIP (gross)" & plotDat1$model == "First-Diff"], 2, plotDat1$lower95[plotDat1$var == "CIP (gross)" & plotDat1$model == "First-Diff"], lty = 2)
segments(2, plotDat1$estimate[plotDat1$var == "CIP (gross)" & plotDat1$model == "First-Diff"], 2, plotDat1$upper95[plotDat1$var == "CIP (gross)" & plotDat1$model == "First-Diff"], lty = 2)
segments(1.5, plotDat1$estimate[plotDat1$var == "CIP (gross)" & plotDat1$model == "LDV"], 1.5, plotDat1$lower95[plotDat1$var == "CIP (gross)" & plotDat1$model == "LDV"], lty = 2)
segments(1.5, plotDat1$estimate[plotDat1$var == "CIP (gross)" & plotDat1$model == "LDV"], 1.5, plotDat1$upper95[plotDat1$var == "CIP (gross)" & plotDat1$model == "LDV"], lty = 2)
segments(1, plotDat1$estimate[plotDat1$var == "CIP (gross)" & plotDat1$model == "Fixed-Effects"], 1, plotDat1$lower90[plotDat1$var == "CIP (gross)" & plotDat1$model == "Fixed-Effects"])
segments(1, plotDat1$estimate[plotDat1$var == "CIP (gross)" & plotDat1$model == "Fixed-Effects"], 1, plotDat1$upper90[plotDat1$var == "CIP (gross)" & plotDat1$model == "Fixed-Effects"])
segments(2, plotDat1$estimate[plotDat1$var == "CIP (gross)" & plotDat1$model == "First-Diff"], 2, plotDat1$lower90[plotDat1$var == "CIP (gross)" & plotDat1$model == "First-Diff"])
segments(2, plotDat1$estimate[plotDat1$var == "CIP (gross)" & plotDat1$model == "First-Diff"], 2, plotDat1$upper90[plotDat1$var == "CIP (gross)" & plotDat1$model == "First-Diff"])
segments(1.5, plotDat1$estimate[plotDat1$var == "CIP (gross)" & plotDat1$model == "LDV"], 1.5, plotDat1$lower90[plotDat1$var == "CIP (gross)" & plotDat1$model == "LDV"])
segments(1.5, plotDat1$estimate[plotDat1$var == "CIP (gross)" & plotDat1$model == "LDV"], 1.5, plotDat1$upper90[plotDat1$var == "CIP (gross)" & plotDat1$model == "LDV"])
text(x=c(1,2,1.5),  par("usr")[3], 
     labels = c("Fixed-Effects","First-Difference ","LDV"),  pos = 1, xpd = TRUE, cex = 1.5)
text(y=c(-6,-4,-2,0,2,4,6),  par("usr")[1],
     labels = c(-6,-4,-2,0,2,4,6), adj = .025,  cex = 1.5)
lines(x = c(0.975,1,1.1,1.2,2.05), y = rep(0,5) , lty = 2, col = "red")


par(mar = c(4,2,4,2) + 0.1, mai = c(1,1,.5,.5), oma = c(1,1,1,1))
plot(plotDat1$estimate[plotDat1$var == "Polls"] ~ c(1.5,1,2), ylim = c(-4, 4), pch = 19, cex = 1.25, col = "black", xaxt='n', yaxt = 'n',  xlab = "",
     ylab = "Marginal Effect on EPS" , cex.lab = 1.7, frame = FALSE, xlim = c(0.95,2))
segments(1, plotDat1$estimate[plotDat1$var == "Polls" & plotDat1$model == "Fixed-Effects"], 1, plotDat1$lower95[plotDat1$var == "Polls" & plotDat1$model == "Fixed-Effects"], lty = 2)
segments(1, plotDat1$estimate[plotDat1$var == "Polls" & plotDat1$model == "Fixed-Effects"], 1, plotDat1$upper95[plotDat1$var == "Polls" & plotDat1$model == "Fixed-Effects"], lty = 2)
segments(2, plotDat1$estimate[plotDat1$var == "Polls" & plotDat1$model == "First-Diff"], 2, plotDat1$lower95[plotDat1$var == "Polls" & plotDat1$model == "First-Diff"], lty = 2)
segments(2, plotDat1$estimate[plotDat1$var == "Polls" & plotDat1$model == "First-Diff"], 2, plotDat1$upper95[plotDat1$var == "Polls" & plotDat1$model == "First-Diff"], lty = 2)
segments(1.5, plotDat1$estimate[plotDat1$var == "Polls" & plotDat1$model == "LDV"], 1.5, plotDat1$lower95[plotDat1$var == "Polls" & plotDat1$model == "LDV"], lty = 2)
segments(1.5, plotDat1$estimate[plotDat1$var == "Polls" & plotDat1$model == "LDV"], 1.5, plotDat1$upper95[plotDat1$var == "Polls" & plotDat1$model == "LDV"], lty = 2)
segments(1, plotDat1$estimate[plotDat1$var == "Polls" & plotDat1$model == "Fixed-Effects"], 1, plotDat1$lower90[plotDat1$var == "Polls" & plotDat1$model == "Fixed-Effects"])
segments(1, plotDat1$estimate[plotDat1$var == "Polls" & plotDat1$model == "Fixed-Effects"], 1, plotDat1$upper90[plotDat1$var == "Polls" & plotDat1$model == "Fixed-Effects"])
segments(2, plotDat1$estimate[plotDat1$var == "Polls" & plotDat1$model == "First-Diff"], 2, plotDat1$lower90[plotDat1$var == "Polls" & plotDat1$model == "First-Diff"])
segments(2, plotDat1$estimate[plotDat1$var == "Polls" & plotDat1$model == "First-Diff"], 2, plotDat1$upper90[plotDat1$var == "Polls" & plotDat1$model == "First-Diff"])
segments(1.5, plotDat1$estimate[plotDat1$var == "Polls" & plotDat1$model == "LDV"], 1.5, plotDat1$lower90[plotDat1$var == "Polls" & plotDat1$model == "LDV"])
segments(1.5, plotDat1$estimate[plotDat1$var == "Polls" & plotDat1$model == "LDV"], 1.5, plotDat1$upper90[plotDat1$var == "Polls" & plotDat1$model == "LDV"])
text(x=c(1,2,1.5),  par("usr")[3], 
     labels = c("Fixed-Effects","First-Difference ","LDV"),  pos = 1, xpd = TRUE, cex = 1.5)
text(y=c(-6,-4,-2,0,2,4,6),  par("usr")[1],
     labels = c(-6,-4,-2,0,2,4,6), adj = .025,  cex = 1.5)
lines(x = c(0.975,1,1.1,1.2,2.05), y = rep(0,5) , lty = 2, col = "red")


##################################################################
# Figure 4, panel A first, then panel B                          #
##################################################################

plotDat3 <- data.frame(rbind(
  cbind(model = "LDV", estimate = mod.3.c$coefficients[length(mod.3.c$coefficients)-1], se = mod.3.c.se[10], var = "PM's CIP net of Green Part(ies)"),
  cbind(model = "Fixed-Effects", estimate = mod.3.a$coefficients[length(mod.3.a$coefficients)], se = mod.3.a.se[length(mod.3.a.se)], var = "PM's CIP net of Green Part(ies)"),
  cbind(model = "First-Diff", estimate = mod.3.b$coefficients[length(mod.3.b$coefficients)], se = mod.3.b.se[length(mod.3.b.se)], var = "PM's CIP net of Green Part(ies)")))


rownames(plotDat3) <- 1:nrow(plotDat3)
plotDat3$estimate <- as.numeric(as.character(plotDat3$estimate))
plotDat3$se <- as.numeric(as.character(plotDat3$se))
plotDat3$upper95 <- plotDat3$estimate + plotDat3$se * qnorm(0.975)
plotDat3$lower95 <- plotDat3$estimate - plotDat3$se * qnorm(0.975)
plotDat3$upper90 <- plotDat3$estimate + plotDat3$se * qnorm(0.95)
plotDat3$lower90 <- plotDat3$estimate - plotDat3$se * qnorm(0.95)
plotDat3$var <- as.character(plotDat3$var)
plotDat3$model <- as.character(plotDat3$model)


par(mar = c(4,2,4,2) + 0.1, mai = c(1,1,.5,.5), oma = c(1,1,1,1))
plot(plotDat3$estimate ~ c(1.5,1,2), ylim = c(-1, .5), pch = 19, cex = 1.25, col = "black", xaxt='n', yaxt = 'n',  xlab = "",
     ylab = "Marginal Effect on EPS" , cex.lab = 1.7, frame = FALSE, xlim = c(0.95,2))
segments(1, plotDat3$estimate[plotDat3$model == "Fixed-Effects"], 1, plotDat3$lower95[plotDat3$model == "Fixed-Effects"], lty = 2)
segments(1, plotDat3$estimate[plotDat3$model == "Fixed-Effects"], 1, plotDat3$upper95[plotDat3$model == "Fixed-Effects"], lty = 2)
segments(1, plotDat3$estimate[plotDat3$model == "Fixed-Effects"], 1, plotDat3$lower90[plotDat3$model == "Fixed-Effects"])
segments(1, plotDat3$estimate[plotDat3$model == "Fixed-Effects"], 1, plotDat3$upper90[plotDat3$model == "Fixed-Effects"])
segments(2, plotDat3$estimate[plotDat3$model == "First-Diff"], 2, plotDat3$lower95[plotDat3$model == "First-Diff"], lty = 2)
segments(2, plotDat3$estimate[plotDat3$model == "First-Diff"], 2, plotDat3$upper95[plotDat3$model == "First-Diff"], lty = 2)
segments(2, plotDat3$estimate[plotDat3$model == "First-Diff"], 2, plotDat3$lower90[plotDat3$model == "First-Diff"])
segments(2, plotDat3$estimate[plotDat3$model == "First-Diff"], 2, plotDat3$upper90[plotDat3$model == "First-Diff"])
segments(1.5, plotDat3$estimate[plotDat3$model == "LDV"], 1.5, plotDat3$lower95[plotDat3$model == "LDV"], lty = 2)
segments(1.5, plotDat3$estimate[plotDat3$model == "LDV"], 1.5, plotDat3$upper95[plotDat3$model == "LDV"], lty = 2)
segments(1.5, plotDat3$estimate[plotDat3$model == "LDV"], 1.5, plotDat3$lower90[plotDat3$model == "LDV"])
segments(1.5, plotDat3$estimate[plotDat3$model == "LDV"], 1.5, plotDat3$upper90[plotDat3$model == "LDV"])
text(x=c(1,2,1.5),  par("usr")[3], 
     labels = c("Fixed-Effects","First-Difference ","LDV"),  pos = 1, xpd = TRUE, cex = 1.5)
text(y=c(-1,0,1),  par("usr")[1],
     labels = c(-1,0,1), adj = .025,  cex = 1.5)
lines(x = c(0.975,1,1.1,1.2,2.05), y = rep(0,5) , lty = 2, col = "red")



set.seed(3)
b.sim <- mvrnorm(2000, coef(mod.3.e), vcovHC(mod.3.e,type="HC0",cluster="group"))
dim(b.sim)
colnames(b.sim)

cond.var.scale = seq(min(mod.3.e[["model"]]$pm_cip_eco), max(mod.3.e[["model"]]$pm_cip_eco), length.out = 175)


pred.d0 <- cbind(1, as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "denmark"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "finland"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "germany"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "ireland"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "italy"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "netherlands"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "portugal"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "sweden"),
                 mod.3.e[["model"]]$cab_minority,mod.3.e[["model"]]$kyoto,
                 mod.3.e[["model"]]$gdp_quarter, mod.3.e[["model"]]$environ_cab_mean,
                 mod.3.e[["model"]]$total_ghg,mod.3.e[["model"]]$eco_gov,
                 mod.3.e[["model"]]$environ_eco, cond.var.scale, 0, mod.3.e[["model"]]$year,  0 * cond.var.scale )


pred.d1 <- cbind(1, as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "denmark"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "finland"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "germany"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "ireland"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "italy"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "netherlands"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "portugal"),
                 as.numeric(mod.3.e[["model"]]$`as.factor(country_id)` == "sweden"),
                 mod.3.e[["model"]]$cab_minority,mod.3.e[["model"]]$kyoto,
                 mod.3.e[["model"]]$gdp_quarter, mod.3.e[["model"]]$environ_cab_mean,
                 mod.3.e[["model"]]$total_ghg,mod.3.e[["model"]]$eco_gov,
                 mod.3.e[["model"]]$environ_eco, cond.var.scale, 1, mod.3.e[["model"]]$year, 1 * cond.var.scale )            

preds.0 <- pred.d0 %*% t(b.sim)
preds.1 <- pred.d1 %*% t(b.sim)

mfx.d.10 <- apply(preds.1 - preds.0, 1, 
                  quantile, c(.5, .025, .975))

par(mar = c(4,2,4,2) + 0.1, mai = c(1,1,.5,.5), oma = c(1,1,1,1))
plot(cond.var.scale, mfx.d.10[1, ], xlim = c(0,1), ylim = c(-4,7), type = "n", ylab = "AME of Green Party's gross CIP",
     xlab = "PM's net CIP excl. Green Party", cex.axis = 2, cex.lab = 2, frame = F, xaxt = "n", yaxt = "n") 
polygon(c(cond.var.scale, rev(cond.var.scale)), y = c(mfx.d.10[2, ], rev(mfx.d.10[3, ])),
        col = 'gray80', border = NA)
lines(cond.var.scale, mfx.d.10[1, ]);lines(cond.var.scale, rep(0,length(cond.var.scale)), col = "red", lty = 2)
rug(mod.3.e[["model"]]$pm_cip_eco, side = 1)
text(x = c(0,0.2,0.4,0.6,0.8,1), par("usr")[3], labels = c(0,0.2,0.4,0.6,0.8,1),  pos = 1, xpd = TRUE, cex = 2)
text(y=c(-4,-2,0,2,4,6),  par("usr")[1],
     labels = c(-4,-2,0,2,4,6), adj = .025,  cex = 2)

